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Abstract 

Averaged  equations  of  motion  were  numerically 
integrated  over  a  time  span  of  ten  Earth  years  for  various 
inclinations,  eccentricities,  and  perigee  heights.  Each 
orbit  (i.e.  set  of  initial  conditions)  produced  a  standard 
deviation  of  the  variations  in  eccentricity  (SDE)  and 
inclination  (SDI) .  These  were  calculated  using  a  polynomial 
approximation  to  the  variations.  Surface  plots  of  SDI  &  SDE 
vs  the  initial  conditions  are  then  examined  to  ascertain  the 
critical  inclinations. 
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NUMERICAL  ANALYSIS  OP  CRITICAL  INCLINATIONS  ABOUT  THE  PLANET 


MARS 

I.  Introduction 

Satellites  are  placed  in  a  variety  of  orbits,  usually 
determined  by  the  mission  requirements  of  the  satellite 
payload.  For  instance,  typical  planetary  mapping  missions 
require  low,  altitude,  high  inclination,  and  nearly  circular 
orbits.  Conversely,  communication  relay  satellites  usually 
occupy  extremely  high  altitude  and  near  zero  inclination 
orbits.  One  such  satellite  is  the  commercial  Intelsat 
satellite  in  Earth  geosynchronous  orbit.  The  basic 
Kepler ian  two  body  equations  are  used  to  find  an  orbit  that 
meets  the  general  requirements  such  as  coverage,  refresh 
rate,  etc.  However,  various  gravitational,  atmospheric,  and 
other  effects  tend  to  perturb  the  orbit  from  its  initial 
conditions.  These  effects  can  be  relatively  large  and 
primarily  affect  the  inclination  and  eccentricity  of  an 
orbit . 

The  amplitudes  of  these  variations  are  primarily 
dependent  on  the  orbit's  initial  eccentricity  and 
inclination.  Higher  initial  eccentricities  show  a  direct 
correlation  to  larger  variations.  The  effect  of  initial 
inclination  on  the  variations  is  more  complex.  Those 
initial  inclinations  that  are  determined  to  cause  large 
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variations  in  eccentricity  and  inclination  are  termed 
critical  inclinations.  Figures  1  &  2  illustrate  the 
variations  in  eccentricity  and  inclination  for  an  example 
orbit. 


B  For  the  majority  of  satellite  missions,  such  as  Earth 

g  communication  satellites,  it  is  desirable  to  minimize  the 

perturbations.  The  initial  orbits  are  selected  to  meet 
^  mission  requirements  and  the  satellites  are  station  kept  in 

those  orbits  by  thrusters.  These  thrusters  require  fuel, 

I  which  is  often  the  limiting  factor  of  satellite  lifetime. 

■  An  alternate  approach  seeks  to  select  orbit  initial 
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conditions  to  maximize  the  perturbations.  Such  an  orbit 
could  be  used  as  a  transfer  orbit  to  save  maneuvering  fuel. 
However,  such  a  transfer  often  takes  much  longer  than 
mission  constraints  allow. 


Whether  the  designer  seeks  to  minimize  or  maximize  the 
perturbations,  the  locations  of  the  critical  inclinations 
are  required.  Therefore,  the  evolution  of  the  selected 
orbit  over  the  lifetime  of  the  satellite  must  be  examined. 
Both  an  analytic  and  numerical  approach  were  applied  in 
finding  the  critical  inclinations. 

Both  approaches  start  with  the  same  governing  equations 
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of  motion,  the  Lagrange  Planetary  Equations  in  their 
disturbing  function  form.  The  disturbing  function  accounted 
for  only  the  predominate  gravitational  perturbations  of  the 
Mar's  equatorial  bulge  (J2)  and  solar  third  body 
perturbations.  The  Mars  J2  zonal  harmonic  is  two  orders  of 
magnitude  larger  than  the  other  harmonics  and  is  twice  as 
large  as  the  Earths.  This  simplified  the  equations  of 
motion  to  permit  an  analytical  analysis  and  reduced  the 
integration  time  for  the  numeric  approach.  The  analytic 
analysis  represents  the  limiting  case  of  variations  over  an 
infinite  time  span.  Therefore,  a  full  numerical  analysis 
was  required  to  find  intermediate  critical  inclinations. 

The  numerical  approach  integrated  orbits  for  six 
different  perigee  heights,  over  a  full  range  of  initial 
eccentricities  and  inclinations.  This  involved  the 
integration  of  over  56,000  orbits  for  a  ten  Earth  year  time 
span.  The  integration  of  such  a  large  number  of  orbits  was 
necessary  to  provide  sufficient  information  for  the 
determination  of  the  critical  inclinations.  However,  the 
shear  amount  of  data  precludes  analysis  of  each  individual 
orbit.  Therefore,  an  approach  involving  the  determination 
of  a  figure  of  merit  for  each  orbit  was  used  to  reduce  the 
amount  of  data  to  a  manageable  size.  The  standard 
deviations  in  eccentricity  and  inclination  were  selected  as 
the  figures  of  merit  and  calculated  for  each  orbit.  Three- 
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dimensional  surface  plots  of  these  deviations  vs  the  initial 
eccentricities  and  inclinations  were  then  used  to  determine 
the  critical  inclinations.  The  affect  of  using  linear, 
quadratic,  and  cubic  data  fits  in  calculating  the  standard 
deviations  was  also  examined.  Figures  1  &  2  show  a  cubic 
fit  to  the  data. 
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II.  Analytical  Development 


Reference  Frame 


A  Mars  centered  reference  frame  and  classical  elements  were 
used  describe  the  motion  of  the  satellite. 


Satellite  orbital  elements  are  defined  as  follows: 
r  =  satellite  radius  vector 
i  =  inclination  of  satellite  orbit 

fl  =  argument  of  ascending  node  for  satellite  orbit 
a)  =  argument  of  periapsis  for  satellite  orbit 
0  =  satellite  central  angle 

By  referencing  everything  to  a  Mars  centered  coordinate 
frame,  the  Sun  appears  to  orbit  Mars  in  this  reference  frame. 
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Therefore,  the  motion  of  the  Sun  with  respect  to  Mars  can  be 
described  by  a  set  of  orbital  elements  as  if  the  Sun  was  orbiting 
Mars.  The  Sun's  pseudo  orbital  elements  are  defined  as  follows: 
rs  =  Sun  radius  vector 
0s  =  Sun  central  angle 
I  =  inclination  of  Sun  orbit 
A  =  right  ascension  of  the  Sun 
L  =  Planet  centered  latitude  of  the  Sun 
Equations  of  Motion 

For  a  basic  two  body  problem  of  a  satellite  orbiting  Mars, 
the  motion  can  be  described  by  the  six  classical  orbital  elements 
that  are  constant  in  the  inertial  frame.  Once  these  are  known 
for  some  epoch  time,  the  position  and  velocity  of  the  satellite 
can  be  determined  for  any  other  time  using  the  two  body 
equations.  In  reality,  the  two  body  elements  are  not  constant. 
They  change  with  time  due  to  the  perturbative  effects  of  drag. 
Mars  geopotential,  Sun  third  body  effects,  etc.  Since  these 
changes  are  relatively  slow,  the  Keplerian  elements  can  be  used 
to  describe  the  satellite  motion  when  combined  with  equations  of 
variation.  The  equations  of  variation  describe  the  change  of  the 
keplerian  elements  with  time,  based  upon  the  perturbations.  This 
method  is  known  as  variation  of  parameters.  The  Lagrange 
disturbing  function  approach  was  used  for  the  equations  of 
variation . 

In  the  Lagrange  approach,  the  acceleration  of  the  satellite 
is  written  as  the  gradient  of  the  gravitational  potentials. 
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Note,  this  is  only  valid  for  conservative  forces  such  as  gravity. 


-^4  =  V(Rm  +  Rs)  (1) 

dt2 

V  =  gradient  operator 

R^  =  Mars  gravity  potential 

Rs  =  Sun  gravity  potential 

Through  Canonical  transformations,  the  two  body  motion  may 
be  eliminated  and  the  Lagrange  equations  for  the  perturbing 
motion  formed.  The  Lagrange  Planetary  Equations  in  their 
disturbing  function  form  (4,476-483): 

e  =  eccentricity  of  satellite  orbit 
i  =  inclination  of  satellite  orbit 
b  =  semi-minor  axis  of  satellite  orbit 
M  =  Mean  Anomaly  of  satellite  orbit 
n  =  mean  motion  of  satellite 
M  =  gravitational  parameter 

(4 . 2828287E04  Kg3/sec2  for  Mars) 

(1. 3271544E11  Kg3/sec2  for  Sun) 

R  =  the  disturbing  function 
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de 

dt 


b2  dR  _  b  djg 

na*e  dM  na2e  d w 

di  _  _  1 _ dR_  +  cos  i  dR 

dt  nab  sin i  dQ  nab  sini  Sco 

da  _  2  dR 

dt  na  dM 

dfl  _  1  aj? 

dt  nab  sini  di 

dco  _  _  cosi  di?  b  dR 

dt  nab  sin i  di  naze  de 

dW  =  n  -  2  dR  -  b2  dR 
dt  na  da  na*e  de 


(2) 


Only  the  Mars  J2  and  Sun  third  body  disturbing  potentials 
were  considered.  Their  development  followed  that  of  Breakwell's 
Investigation  of  High  Eccentricity  Orbits  About  Mars  (5)  and 
Durand's  June  1989  MS  Thesis  (7). 
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Mars  Disturbing  Function 


The  first  non-zero  zonal  harmonic  is  the  J2  term  that 
accounts  for  the  equatorial  bulge.  This  is  the  predominate 
geopotential  term  for  most  planets  and  the  only  harmonic  examined 
in  both  the  analytic  and  numerical  developments.  The  J2 
disturbing  function  has  been  shown  to  be 

r  =  (1  _  3  sin25)  (3) 

*  2r 3 


where  Re  is  the  Mars  equatorial  radius  and  <5  is  the  declination. 

A  transformation  of  the  declination  to  a  form  containing  the  true 
anomaly  and  the  argument  of  periapsis  was  accomplished  using 
spherical  geometry.  From  the  spherical  geometry 

sin5  =  sin i  since 

=  sini  sin  ( iz  -  (  f  +  w  )  )  (4) 

=  sin  i  sin  (  £  +  to  ) 


Squaring  this  equation  yields 


sin2  6  =  sin 2  i  sin 2  (  f  +  o) 

=  —  sin2  ifl  -  cos2(f  +  o))l 
2  1  J 


(5) 
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Figure  4.  Spherical  Trigonometry  of  an  Orbit 


Figure  5.  Geometry  of  the  Orbital  Plane 
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The  transformed  disturbing  function  is  thus 


R 


m 


m^2 

2r 3 


3 

2 


si  mi  + 


sin2  i  cos  2  (  f  +  o> ) 


(6) 


Averaging  of  the  Mars  Disturbing  Function 

The  mean  anomaly  was  averaged  out  of  the  Mars  disturbing 
function  by  integrating  the  function  over  one  orbit  and  dividing 
by  the  change  in  mean  anomaly  (27 r) . 

r ^  (?) 

2  7t  Jo 


The  integration  was  accomplished  using  Hansen's  Coefficients  (2). 
A  brief  description  and  derivation  of  this  is  included  in 
Appendix  A.  Substituting  in  the  disturbing  function  yields 


R, 


l  r2lt  H  m  Rp  i  3  .  ,  , 

=  —  /  2  6  1  -  —  sin2  i 

2n  Jo  2r 3  \  2 


dM 


(8) 


J-f 

2u  Jo 


271  3  M>  ^ 


-J2— 2 — 2.  sin2  i  cos  2  (  f+o> )  dM 
4r3 


Let 


A  = 


j  1 

2a3  \ 


dM 


B  = 


4a3 


;in2i  {  — —  f2n  (—)  3  cos2(f+co)  dM 
(  2n  Jo  \  a  / 


(9) 
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Using  trigonometry,  B  was  transformed  to 


S  - 


B  = 


3 


4a3 
3  P>  2 


4a- 


sin2i  1-i-  f2n  (—]  (cos2f  cos2w  -  sin2f  sin2o> 
(  2  m  Jo  \  a  / 

sirri  |  (cos2o>)  — f2*  f— \  3  cos2 f  dM  + 

1  2m  Jo  \  a/ 


dM 


i  /*2it  /  r  \-3 

sin2o>  -  /  —  sin2f  aM 

2m  Jo  \  a  / 


(10) 


By  employing  Hansen's  coefficients  from  Appendix  A 


1 

pin 

I  —  )'3  dM  =  Xo3,0  =  (l-e2)'3/2 

2m 

Jo 

\  a  / 

1 

f  2n 

(  — )  3  sin2f  dM  =  0 

2  m 

Jo 

\  a  / 

1 

r2n 

(  —  T3  cos  2  f  dM  =  X0'3’2  =  0 

2m 

Jo 

\  a  / 

From  this  "B"  is  found  to  equal  zero.  Therefore  the  averaged 
disturbing  function  is  equal  to  the  "A”  equation. 


477 


M1  m  ^2  Re 


.  (l  -  —  sin2  i 
2a3  (1-e2) 3/2  \  2 


(12) 
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Solar  Disturbing  Function 

The  third  body  effects  of  the  Sun  on  the  satellite  motion 
may  be  derived  from  Newton's  basic  Law  of  Universal  Gravitation. 


F  = 

-  G/m 

m2 

r*  + 

^  ~ 

rl 

rL 

X21  + 

r3  31 

X31 

F  - 

-G/n2 

'  mi 

F12  + 

%  —  | 

”  2 

3  ^2 
r32 

<t2rx  _ 

A  = 

-  c 

m0  _ 

z  -r 

dt2 

wi 

rh 

21  3 

t3l 

d2F,  _ 

II  = 

-  G 

m0  - 

-  +  ™3 

dt2 

m2 

r?/ 

12  r3 

r32 

■31 


■  32 


(13) 


(14) 


d2r 


21 


d2r ;  d2Fz 


d  t‘ 


dt: 


dt‘ 


d  t~ 


G 


(m. 


z-. 


-r  m.)  _ 

“  ~  rzi 


Gmz 


21 


•  32 


31 


r33i 


d2T,l 

-  -  (Z 

^  J  + 

^  —  1 

+  G 

F  + 

m3 

dt2 

3  21 

X21 

r3  31 

X31 

3  ^ 

L  Zl2 

rir"J 

(15) 
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Since  m2  >>  ml,  let  m2  +  ml  =  m2 


d2r 
u  -‘-21 

dt2 


(16) 


With  a  change  of  variables  to  the  notation  of  figure  3. 


The  acceleration  of  the  satellite  is 


d2r 

dt2 


V(Rm  +  Rs) 


(18) 


The  first  term  of  the  above  equation  is  the  Mars  and  satellite 
two  body  problem.  The  second  term  is  the  acceleration  due  to  the 
Sun  third  body  effects.  As  shown  above,  the  acceleration  of  the 
satellite  can  be  expressed  as  the  gradient  of  the  gravitational 
potentials  of  the  Sun  and  Mars.  Therefore,  the  gradient  of  the 
Sun's  potential  is 


VRS  =  n 


s 


(19) 


By  finding  the  equation  that  satisfies  this  condition,  the  Sun's 
gravitational  potential  is  found  to  be 
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R  < 


r-r. 


rs  ~  r 


(20) 


Now  this  must  be  put  in  the  form  of  a  disturbing  function  to  the 
Mars-satellite  two  body  problem. 


P  =  rs  -  r 

P2  =  \rs  -  r\2  =  (rs  -  F)  •  (rs  -  r) 
P2  =  rs*Fs  -  2r-rs  +  t-r 
p2  =  r2s  -  2 rrs  cos B  +  r2 
r^r” 

where  B  =  - - 


(21) 


By  factoring  out  rs2  and  taking  the  square  root 


P  =  r. 


1  + 


2r 


cos  fl 


1/2 


=  r, 


(22) 


Substituting  in,  the  Mars  gravitational  potential  becomes 


Rs  = 


1  +  - 


2r 


\-l/2 


cos  B 


rcosB 


(23) 
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Using  the  binomial  expansion 


1  + 


2  r 


r 2 
\  Zs 


COS  B 


-1/2 


/]  =  0 


.1  ^ 
2 


n 


2  r 


COS  B 


/ 


(24) 


1 

—  -  —  cos  B 

+  2 

'^-i£coss) 

2 

r|  r« 

8 

[4  j 

The  higher  order  terms  above  (r/rs)2  can  be  neglected  since 
rs  >>>  r.  The  above  series  is  approximated  as 


1  + 

f-£i-i£coss] 

1-1/2 

=  i-2 

'll-l£cosB] 

3 

+  — 

4r2COS2B) 

2 

U  J 

8 

W  J 

The  Sun  gravitational  potential  becomes 


= 


I c 


1  - 


1  + 


r  *  -Ml  cos 2B 
2 


2rJ 


2  r] 


2  r\ 


(  3  cos 2 B  -  1 ) 


(26) 


The  Suns  gravitational  potential  at  Mars  is  equal  to  the  two  body 
term  plus  the  disturbing  function. 


RS  = 


+  R* 


(27) 
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Therefore,  the  solar  disturbing  function  is 


"  3 

2  r3s 

[3  COS2  3 

- 1] 

\xsr2 

2 

|  -  T 

2  r\ 

[  1  rrs  J 

u_ 

The  radial  vectors  are  now  transformed  to  the  orbital 
elements  of  the  satellite  and  Sun,  defined  in  figure  (1) .  For 
ease  of  notation,  the  sin  and  cos  functions  will  be  abbreviated 
using  S  =  sin  and  C  =  cos.  From  the  physical  geometry,  the  unit 
vectors  in  the  radial  satellite  and  Sun  directions  are 


''■7 


= 


CqC0  Sq  Ci  Sq 
SqCq  +  CqC2S8 

SiS0 


C\CL 


(29) 


The  dot  product  of  these  unit  vectors  is 


er'es  = 


r  *r£ 
rr„ 


CqC9CaCl  -  Sa C i Sq CA C,  +  +  CQCi£'eCL5A  +  S L 

~  ^  Cq 0^  +  SqSa)  +  Ci56Ct(CQ5A  -  SqCa)  +  S^SqS^ 

=  cqclc(Q  _  A)  +  CiS6CLSl Q  .  A)  +  SiSeSL 


(30) 
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Substituting  into  the  solar  disturbing  function  results  in 


\isr- 


Rs  - 


2  r; 


{3  [  CeCLC(Q  .  A)  +  CiSeCLS{ q  .  a)  +  2  i}  (31) 


Expanding  squared  term  yields 


*5  =  -^“M^Cg2  +  2  K2  Cq Sq  +  ^S2] 
2r2 


1} 


(32) 


Where 


*i  = 


= 


= 


(  CLC(  q  -  a  )  )  2 

SLCLSi  C(Q  -  A)  "  ^i5(Q  -  A)  ^(Q  -  A) 

<<W?(Q.A)  -  SLSi)2 


(33) 


The  disturbing  function  is  now  written  in  terms  of  the  satellites 
true  anomaly  using 

0  =  f  +  0) 


Ce  =  (1  +  C20) 

Se  =  j  (1  -  C20) 
2  iSg  Cg  =  S2  g 


(34) 


The  final  form  of  the  solar  disturbing  function  is 


= 


^sr  I  3 


2  r: 


|  2  ^2  +  ^  ^2 


+  2^  _  i 

2 (ftu)  2  3  X 


(35) 
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Where 


Ai  =  *1  -  *3 

^2  *  *2  (36) 
a3  =  jq  +  k2 
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Averaging  the  Solar  Disturbing  Function 


Just  as  with  the  Mars  gravitational  disturbing  function,  the 
fast  periodic  mean  anomaly  can  be  averaged  out  of  the  solar 
disturbing  function.  This  is  accomplished  by  integrating  the 
disturbing  function  with  respect  to  the  mean  anomaly  over  one 
orbit  and  dividing  by  the  total  change  in  mean  anomaly  (27 r) . 

Thus  the  average  solar  disturbing  function  is 

F  =  -L  [2*Rs  dM  (37) 

s  2k  h  s 


Before  integrating,  the  true  anomaly  must  be  separated  from  the 
argument  of  periapsis  o .  Using  the  trigonometric  relations 


"  !~2f('2<o  S2fS2(ji 


~  FfFw  +  F/^2  <j 


(38) 


The  solar  disturbing  function  may  now  be  written  as 


= 


(  3 


2  r 


3  2 


-  ^iC2{C2w  2AxS2fS2u>  +  3A2S2fC2a)  + 


+  ~$A2C2fS2U)  *■  ^  A3  i 


(39) 


f /  3 


2r'  n  * 


.  3 A2s2^c2f  ♦  (3A2C2u  -  |^(J)S2f  * 


*  ~A,  -  1 
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Upon  rearranging  to  a  form  useful  for  employing  Hansen's 
Coefficients 


■  c- 

s-  *  (-h  -  'Ut  I 


Employing  Hansen's  coefficients  (Appendix  A)  to  average  out  the 
Mean  Anomaly  using 


£  r  (if  ^  - 

1  r2n  I  r\2 

— —  /  —)  S2f  dM  =  0 

2  k  Jo  \  <a  / 


2,2  _  5e2 
0  - 


After  substitution,  the  averaged  solar  disturbing  function 
becomes 


ItAn  ^  a.  *»  a  o  \  5e^  +  I  —  A-,  -  l\(l  +  —  e2 


3  2 


- Ai  C2(J  +  3 A2S2(i 


j^{JTe2AlC‘-  +  15e2^52u  -  (2  +  3e2)(|A3  -  1 
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Substituting  in  the  "A"  variables  results  in 


Rs  =  ~e2[c£2C?0_A)  -  -  S^)2]^ 

4  rs  l 

+  15  e2  |CLSi5i  C(Q.A)  -  CL  Ct  C(Q.A)  *?(Q_A)  }  S2(i)  + 

+  (2  +  3e2)[|[C'C<Q-A)  +  (CiC^(Q- A,  -  Sl^)2]  -  l]  } 


However,  this  form  still  contains  periodic  effects  due  to 
the  motion  of  the  Sun.  These  may  be  averaged  out  when  the  mean 
motion  of  the  Sun  is  greater  than  the  rate  of  change  of  either  w 
or  Cl.  To  find  the  condition  under  which  this  holds,  use  the 
Lagrange  Planetary  Equations  to  get  the  changes  in  w  and  ft.  Use 
only  the  J2  portion  of  the  disturbing  function  (Rm)  to  simplify 
the  comparison.  This  is  valid  since,  the  affect  of  the  J2 
disturbing  function  is  much  greater  than  the  solar  effects.  The 
necessary  partial  of  the  Mars  disturbing  function  are 


3*. 

di 


de 


ln2J2Rl 
2  (1-e2) 3/2 


sini  cos i 


3  n 2  J2  Rl  e  i 
2  (1-e2) 5/2  \ 


—  sin 


2 


(44) 
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The  rates  of  change  of  n  and  o  are  then 


do) 

dt 


dQ  =  _  3  nJ2Rl 

dt  2a2  ( 1-e2) 2 


COS  1 


3  nJ2Rg 


2  a2  (1-e2) 


2  \  2 


i1 '  isi 


sin2i  +  cos2  i 


(45) 


By  comparing  the  magnitudes  of  these  rates  of  change 

dQ  _  do  w  J2Re  ^m/2 
dt  dt  a7/2  (l-e2)2 


Therefore,  the  averaging  of  the  solar  motion  is  valid  when 


n. 


p2i.1/2 

lKeV-m 


a1’2  (1-e2)2 


(47) 


This  occurs  for  high  eccentricity  orbits  with  the  lower  bound  of 
eccentricity  decreasing  with  increasing  periapsis  radius. 

The  majority  of  the  orbits  analyzed  in  the  numerical  section  fall 
into  this  category. 

After  several  change  of  variables  and  assuming  the  Sun's 
pseudo  orbit  to  be  circular,  the  time  varying  solar  terms  can  be 
removed.  The  slowly  varying  solar  disturbing  function  becomes 
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*s  = 


n;  a* 


[2  +3  e2] 
3 


Tc“s|*sJ  *  \clSlslSlica 


1  -  —Si 

4  2  1 


+  -i^e2 

4 


JA5^(1  +  c|)  (c20+2w+c20_2j 


r,)(3  C2i- 1)} 

-  CzSzCImSIm  <C0t2tt  +  C0_aJ  h-  Sf(l  - 


15 


-  —  e2  \—CTSj  (C, 

O  /I  x  -*s  ^ 


Q-2u  ^2Q*2u 


)  2  ^isSlCI''  (Cq_2(i)  Cs 


Q*2co 


(48) 


where  the  mean  motion  of  the  Sun  is 


The  Combined  Mars  and  Solar  Disturbing  Function 

The  mean  anomaly  and  the  periodic  effects  of  the  Sun's 
motion  have  been  successfully  averaged  out  of  the  disturbing 
functions,  resulting  in 

R  =  + 


R 


J2  K 


2 a3  ( 1-e2) 


2  \  3/2 


1  -  —  sin2  i 
2 


n2s  a2 


[2  +  3  e2] 


-|c2flSj.Sf  -  -|  0,3,5,  c.ca 


*  K1  "  e‘  {-ysf,<i-cf>  (c2q.2»*c2Q_2l 


15 


-  Ci^jCrS.  <CQ.2„*C0.2„)  *  S?(l  -  4s|,lc„ 


7CISIj^2D-2o  ^2  "7  Sx'SjCj-*  (Cq-2o  03-f2o>) 


(50) 


(51) 
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Determination  of  Critical  Inclinations  by  Analytic  Methods 


The  long  period  variations  in  eccentricity  and  inclination 
may  now  be  obtained  chrough  use  of  Hamiltonian  mechanics  and 
canonical  transformations.  Detailed  explanitions  are  in  ref  (5). 

The  transformation  of  the  disturbing  function  results  in  an 
eguation  with  six  divisors.  Eleven  critical  inclinations  result 
when  these  are  equated  to  zero. 


Table  1.  Critical  Inclinations 


63.4  deg 

116.6  deg 

46.4  deg 

106.8  deg 

73.2  deg 

133.6  deg 

56.1  deg 

111.0  deg 

69.0  deg 

123.9  deg 

90.0  deg 

These  inclinations  represent  the  limiting  case  for  critical 
inclinations  and  were  derived  with  the  simplification  of  solar 
averaging.  Due  to  the  complexity  of  the  problem,  a  systematic 
numerical  approach  is  necessary  to  ascertain  the  true  nature  of 
the  critical  inclinations  without  the  above  simplifications. 
However,  the  results  in  table  (1)  should  roughly  correspond  to 
any  numeric  results. 
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III.  Numerical  Development  and  Results 

The  analytical  development  produced  eleven  critical 
inclinations  resulting  from  resonant  situations.  This 
occurred  when  terms  in  the  denominator  of  the  equations  of 
motion  went  to  zero,  thus  indicating  large  amplitude  motion. 
However,  in  reality  these  amplitudes  are  not  infinite.  They 
build  up  over  time  due  to  the  aaaitive  effect  of  the 
resonance  condition.  The  equations  of  motion  are  still 
valid  and  can  be  numerically  integrated  as  long  as  the 
integration  time  is  small  enough  that  the  amplitude  does  not 
exceed  the  limitations  of  the  computer  (overflow  condition) . 
The  effect  of  the  resonance  will  still  be  visible  even  with 
the  limited  integration  time.  The  amplitude  of  the  motion 
will  not  go  to  infinity,  but  will  be  large  compared  to  non¬ 
resonance  conditions.  Therefore,  a  numerical  approach  may 
be  used  to  find  the  critical  inclinations. 

This  involves  integrating  a  sufficient  number  of  orbits 
over  a  range  of  eccentricities,  inclinations,  and  periapsis 
altitudes  (initial  conditions)  to  permit  the  determination 
of  critical  inclinations. 

Equations  of  Motion 

The  numerical  analysis  used  a  transformed  set  of  the 
Lagrange  Planetary  Equations  used  in  the  analytic  analysis. 
The  change  of  variables  was  used  to  simplify  the 
implementation  of  the  equations  on  the  computer.  The  same 
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disturbing  function  used  in  the  analytic  approach  was  used 
for  the  numerical  integration.  As  discussed  previously,  the 
mean  anomaly  is  averaged  out  of  the  disturbing  function 
while  still  preserving  the  effects  of  the  resonance 
conditions.  This  averaging  out  of  the  mean  anomaly  eases 
the  numerical  integration  by  simplifying  the  equations  of 
motion.  The  integration  time  is  also  greatly  reduced.  A 
step  size  of  days  instead  of  minutes  can  be  used. 

Since  the  mean  anomaly  was  averaged  out  and  the  primary 
interest  was  the  slow  variable  changes,  such  as  inclination 
and  eccentricity ,  the  dM/dt  equation  was  not  needed. 
Furthermore,  the  semi-major  axis  remains  constant  (da/dt=0) 
since  the  mean  anomaly  was  averaged  out  of  the  disturbing 
function  (dR/dM=0) .  This  reduced  the  number  of  relevant 
Lagrange  Equations  to  the  following  four: 


de 

b2  dR 

b  dR 

dt 

naie  dM 

na3e 

di 

1  dR  .  COS  2 

dR 

dt 

nab  sini  dQ  nab  sini 

d(i3 

dQ 

1  dR 

dt 

nab  sini  di 

(52) 

d<s> 

cos  i  dR 

b 

dR 

dt 

nab  sini  di 

na3e 

de 

where 

b  —  a\j  l -  e 2  Sc 

n-\ 

JL 

a3 
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A  change  of  variables  from  the  classical  (a,e, i,fi, <*>,M) 


to  the  new  variables  (a ,  h ,  k,  fi ,  ,  i)  was  then  used  to 

facilitate  the  integration. 

h  =  e  sinco 
k  =  e  cos  o> 

(53) 

A*  N  ~  ^  +  ( ••  4> ) 


The  transformation  to  the  h  and  k  variables  avoids  the 
singularities  due  to  the  classical  variable  set  (a,e,i,ft,M). 
The  transformed  equations  are  non-singular  for  zero 
eccentricity,  but  singular  for  zero  inclination.  The 
stroboscopic  mean  node  (A.N)  transformation  permits  the 
averaging  of  the  motion  due  to  the  fast  mean  anomaly 
variable,  while  maintaining  the  longer  term  resonance 
effects  as  discussed  in  the  analytic  section.  The  following 
derivatives  are  required  in  the  transformation: 


dR  _  dR  dh  +  dR  dk  _  i^^R  _  ,  dR 

du )  dh  d(o  dk  d co  dh  dk 

dR  _  dR  dh  +  _dR  dk  _  _h  _dR_  +  k  dR 

de  dh  de  dk  de  e  dh  e  dk 


(54) 
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dh  _  dh  de_  +  dh  doo  _  _h_cfe  +  y  do) 
dt  de  dt  00)  eft  e  dt  eft 


dh 

=  -£(- 

b 

1  (k^ 

-  h^\ 

dt 

•l 

na2e 

1  \  dh 

h  dk) 

+  k 

cos  i 

J*  + 

b  I  h 

dR  (  k 

dR\ 

T  A 

|  nah 

sini 

di 

na2e  \  & 

dh  e 

dk) 

eft?  _  i?  8i?  _  ie  cos  i  8i? 
dt  na3  3ic  nah  sini  8i 


die  _  3ic  de  +  3ie  do)  _  _k  de  _  ^  do> 
dt  0e  dt  0o)  dt  e  dt  dt 


dk  _  k(  b  \  I,  dR  v  0i?\ 
dt  '  ef  ^a3ej  T  ~  ^  0iej 


-  h 


cos  i 

i*  + 

b  1  h 

dR  t  k 

dR  \ 

nab  sini 

di 

na2e  \  e 

dh  e 

dk ) 

dk  _  _  b  dR  cosi  8i? 

dt  na3  dd  nab  sini  0i 


(56) 


The  four  relevant  Lagrange  Planetary  Equations  are  then 
transformed  to  the  new  elements. 


dh  _  b  dR  _  k  cot i  8i? 
dt  naz  dk  nab  di 


di 

cot  i  1 

'  k  dR 

-  h  —  \ 

1 

0i? 

dt 

nab  \ 

,  0h 

h  dk) 

nab  sini 

0Q 

dk 

b 

dR  + 

h  cot  i 

0tf 

dt 

na2 

dh 

nab 

0i 

dQ  _  1  3i? 

dt  nab  sin  i  3i 

where  b  =  a  Jl-e2  &  e  =  \/d2  +  k2 


(57) 
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Now  the  disturbing  function  must  also  be  transformed  to 


the  new  variables  (a,h,  k,n,  i)  . 


3 


Li  J Q  I  .  m 

R  =  - — — — - 1  -  — sin^r 

2a3  (l-h2-k2)3/2\  2 


V-s  a  f  15 


4  r: 


^[C^CfQ.A)  -  (C.C^o.A,  -  5L5i)2](k2  -  h 2)  (S8) 


+  3  0[CiSLSiC(Q.A)  -  C£,CiC(Q.A)5(Q_A)]i2ir 
+  (2  +  3 h2  +  3k2)  [-|[cfC?Q_A)  +  [CiCLSi 


Q-A)  ^ 


i)2]  -  i] } 


The  following  four  partial  derivatives  of  this 
disturbing  function  were  then  used  in  the  transformed 
Lagrange  Planetary  equations. 


dR  _  _ 3 

5k  2a3  {l-h2-k2) 5/2 
[isa2k 


(l  -  -j  sin2  ij 


4r; 


{l5[cL2cfQ_A)  -  (CiC^jo.A,  -  Sl5,)2; 


3  0[CL5L5iC(Q_A)  -  Cz>CiC(Q.A)5(Q.A)]  k 
+  6k[|[c2C?Q_A)  +  (C^S^  -  SLSi ) 2]  -  1 


(59) 
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dR  3  ll  mJ0Reh  i  3*o- 

=  - rj-2— g - i  -  —  sm2  j 

dh  2a3  (l-h2-k2) 5/2 1  2 


\isa2h 


{  -  15  [cL2cfQ.A)  -  (C2Cl5(Q.A)  -  5L5i)2] 


+  30[C£SIS,iC(Q_A)  -  CiCjC^.^SjQ.A)]  ^ 

+  6 h  [|[C2cfQ.A)  +  (C&S^  -  S^,)2}  -  1 


dR  _  ~3^mJ2Re 
di  4a3  ( 1  -h2-k2) 


sin2i 


a' 


{15  (SiCLSi q-A)  +  S ,LCi)  (ic2  A  ) 


+  3 +  C’LSiC(Q_A) 5{q.A)]  hk 

-  3  ( 2  +  3h2  +  3k 2)  (C^S,^  -  SLSi)  q-a)  +  SLCi)  } 


^sa  f  15  r. 


[Q.^2  (Q-A)  +2  (  ~SiPi)  ^i^L^fQ-A)]  ^  ^ 


3  0  CL5L5i5(Q.A)  +  CLCiC2  <n_A)  ] 


—  (2  +  3  h2  +  3  k2)  ClS2( q_a,  2  ( C'iC'£,*?(Q-A)  *~iCLCt q_A) 
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These  partial  derivatives  of  the  disturbing  function 
combined  with  the  modified  Lagrange  Planetary  Equations  were 
then  integrated  for  each  set  of  initial  orbital  elements. 

The  full  set  of  transformed  Lagrange  Planetary 
equations,  without  the  simplifications  made  previously,  are 
as  follows  (17,3-3): 
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da 

2  _JL_ 

aj? 

dt 

na  S0 

axw 

dh 

v/l-e2 

ai? 

h  coti  3i? 

dt 

na2 

dk 

na2  \J\-e2  dl 

di 

cot 

i  { 

k  M  -  h  M  + 

dt 

na2  /I 

- e 2  1 

ah  ah 

h  y/l-e: 

na2  Sr. 

1  dR 
S0  dk  w 


P' 


dR 

dk 


N 


dt 

dCi 

dt 

dk 


na *  dh 
1 


na2  y/l-e2  sin  i 
dR 

na 2  \Jl-e2  ^ 


dR  +  3i? 


ao 
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N 


e  =  JJP  +  k2  P'  = 


1+v/l-e2 

(63) 


XN  =  stroboscopic  mean  node 

e 

=  eccentricity 

0  =  Greenwich  hour  angle 

a 

=  semi-major  axis 

n  =  longitude  of  ascending  node 

M 

=  mean  anomaly 

u  =  argument  of  periapsis 

i 

=  inclination 

S0  =  ratio  that  approximates  the  number 
of  orbits  per  planet  revolution 
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Determination  of  Critical  Inclinations 

Critical  inclinations  were  defined  as  those 
inclinations  that  produce  large  variations  in  one  of  the 
orbital  elements.  This  numerical  approach  was  looking  for 
critical  inclinations  in  both  eccentricity  and  inclination. 
Search  Method 

The  search  was  accomplished  by  giving  each  orbit  two 
figures  of  merit.  These  were  the  standard  deviation  in 
eccentricity  and  inclination.  While  using  the  standard 
deviation  as  a  search  criteria  for  critical  inclinations  was 
not  an  ideal  approach,  it  proved  to  be  a  satisfactory  method 
to  handle  the  huge  quantities  of  data  that  required 
processing. 

Procedure 

The  averaged  equations  of  motion  for  a  satellite  were 
numerically  integrated  over  a  10  Earth  year  period  for  a 
given  set  of  initial  conditions.  The  variable  initial 
conditions  were  systematically  incremented  to  cover  the 
entire  spectrum  of  possible  orbits  about  Mars.  Initial 
conditions  included  the  following: 

rp  =  variable  periapsis  radius  (5000  km,  5500  km, 

6000  km,  6500  km,  7000  km,  7500  km) 
e  =  variable  eccentricity  (0.4  to  0.9,  A  =  0.02) 

i  =  variable  inclination  (0.25  to  90.,  A  =  0.25) 


n 


0 


longitude  of  the  ascending  node 


<i)  =  0  argument  of  periapsis 

M  =  0  mean  anomaly 

For  each  orbit  (set  of  initial  conditions) ,  the 
standard  deviation  of  the  variations  in  eccentricity  and 
inclination  were  calculated  using  a  linear,  quadratic,  and 
cubic  polynomial  approximation  to  the  variations. 

Two  dimensional  plots  of  the  standard  deviation  in 
eccentricity  (SDE)  and  standard  deviation  in  inclination 
( SDI )  vs  inclination  were  generated  for  each  initial 
eccentricity  value.  All  the  SDE  vs  inclination  and  SDI  vs 
inclination  plots  were  combined  into  three  dimensional 
surface  plots.  These  surface  plots  were  then  used  to 
identify  the  critical  inclinations  in  both  eccentricity  and 
inclination. 

Surface  Plots 

A  total  of  six  sets  of  three  dimensional  surface  plots 
were  generated.  Each  set  contains  six  plots  of  the  standard 
deviation  in  eccentricity  or  inclination,  for  a  given  data 
fit,  versus  initial  insertion  orbit  eccentricity  and 
inclination.  Linear,  quadratic,  and  cubic  data  fits  were 
used  in  determining  the  standard  deviations. 

Only  the  5000  km  and  7500  km  periapsis  radius  cases 
include  the  linear,  quadratic,  and  cubic  data  fits  for  SDE 
and  SDI.  These  demonstrate  the  effect  of  using  higher  order 
polynomials  in  computing  the  standard  deviations.  Note  the 
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smoothing  of  the  features  and  reduction  in  magnitudes.  The 
remaining  surface  plots  are  for  a  cubic  data  fit.  The  other 
linear  and  quadratic  data  fit  plots  are  in  Appendix  B  and  C. 

Standard  Deviation  in  Eccentricity  (SDE)  Surface  Plots 
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Figure  6.  Standard  Deviation  in  Eccentricity  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5000  Km,  Linear  Data 
Fit 


39 


Figure  7.  Standard  Deviation  in  Eccentricity  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5000  Km,  Quadratic  Data 
Fit 
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Figure  8.  Standard  Deviation  in  Eccentricity  vs  Inclination 
and  Eccentricity:  Periapre  Radius  =  5000  Km,  Cubic  Data  Fit 
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Figure  9.  Standard  Deviation  in  Eccentricity  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5500  Km,  Cubic  Data  Fit 


I 
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Figure  10.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  6000  Km, 
Cubic  Data  Fit 
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Figure  11.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity :  Periapse  Radius  =  6500  Km, 
Cubic  Data  Fit 


Figure  12.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  7000  Km, 
Cubic  Data  Fit 
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Figure  13.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  7500  Km, 
Linear  Data  Fit 
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Figure  14.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  7500  Km, 
Quadratic  Data  Fit 


Figure  15.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  7500  Km, 
Cubic  Data  Fit 
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Standard  Deviation  in  Inclination  fSDI)  Surface  Plots 
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Figure  16.  Standard  Deviation  in  Inclination  vs  Inclination 


and  Eccentricity:  Periapse  Radius  =  5000  Km,  Linear  Data 


Fit 
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Figure  17.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5000  Km,  Quadratic  Data 
Fit 
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figure  18.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5000  Km,  Cubic  Data  Fit 
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Figure  19.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5500  Km,  Cubic  Data  Fit 
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Figure  20.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  6000  Rm,  Cubic  Data  Fit 


54 


Figure  21.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  6500  Km,  Cubic  Data  Fit 
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Figure  22.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  7000  Km,  Cubic  Data  Fit 
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Fit 
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Figure  24.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  7500  Km,  Quadratic  Data 
Fit 


Figure  25.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  7500  Km,  Cubic  Data  Fit 
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I 
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2-D  Critical  Inclination  Plots 


For  each  of  the  SDE  surface  plots,  a  numerical 
procedure  was  used  to  sort  the  data  and  locate  the 
predominate  local  maximums.  These  form  the  ridges  in  the 
surface  plots  and  indicate  the  critical  inclinations.  Cubic 
polynomials  were  fitted  to  the  local  maximum  data  and  then 
plotted  on  a  two  dimension  graph  of  the  initial  eccentricity 
vs  initial  inclination  (orbit  insertion) .  These  are  the 
plots  of  critical  inclination  as  a  function  of  insertion 
eccentricity  and  inclination.  Two-dimensional  contour  plots 
of  the  surfaces  for  a  cubic  data  fit  are  also  included.  The 
linear  and  quadratic  data  fits  not  included  here  are  in 
Appendix  D.  The  features  in  the  SDI  surface  plots  were  not 
predominate  enough  to  locate  accurately  and  were  not 
plotted. 

The  curves  for  the  linear,  quadratic,  and  cubic  data 
fits  in  calculating  the  standard  deviations  were  roughly 
similar.  As  expected,  some  of  the  curves  disappeared  with 
the  higher  order  fits  due  to  the  decreasing  magnitudes  of 
the  standard  deviations.  However,  the  curves  were  also 
shifted  in  position.  They  exhibited  a  shift  to  lower 
eccentricities  and  a  shift  away  from  a  dividing  line  of  53.3 
degrees  initial  inclination. 

Critical  Inclination  Curves  for  SDE 
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Figure  26.  Critical  Inclination  Curves: 
5000  Km,  Linear  Data  Fit 


Periapse  Radius  = 
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Cubic  Data  Fit 
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Figure  32.  Critical  Inclination  Curves:  Periapse  Radius  = 
6000  Km,  Cubic  Data  Fit 


Figure  33.  2-D  Contour  Plot:  Periapse  Radius  =  6000  Km, 

Cubic  Data  Fit 
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Figure  34.  Critical  Inclination  Curves:  Periapse  Radius  = 
6500  Km,  Cubic  Data  Fit 


Figure  35.  2-D  Contour  Plot:  Periapse  Radius  =  6500  Km, 

Cubic  Data  Fit 
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CRITICAL  INCLINATION  CURVES 

E°  -  ?500  KM 

SDE  FOR  LINEAR  Fit 


Figure  38.  Critical  Inclination  Curves:  Periapse  Radius  = 
7500  Rm,  Linear  Data  Fit 
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Figure  42.  Primary  Critical  Inclination  Curves  for  SDE 
(Cubic  Data  Fit) 


Equations  for  curves  of  figure  42: 

i.  =  -1956.28e3+4150.69e2-3008.36e+794.74 


For  0.68  £  eA  z  0.88 
-1019 . 1C  e3  +  1968 ,60e2-1336 . 05  e+367 .94 
For  0.61  £  eB  £  =  0.86 
-  445 . 74  e3  +  570. 55  e2-  233.95e+  80.66 
For  0.59  •s  ec  £  0.83 
-1121 .95  e3  +  2053 .21  e2-1324 .73  e+344 .94 
For  0.52  £  eD  £  0.79 
-7  51.01e3  +  1235.47  e2-  743 . 30  e+206 .72 
For  0.46  £  eE  <.  0.77 
-731 . 68  e3  +  110  8 . 97  e2  -6 19 . 54  e+17 0 . 27 


For  0.40  £  eF  £  0.74 
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IV.  Conclusions  and  Recommendations 


Critical  inclinations  were  defined  as  those  initial 
inclinations  which  result  in  a  large  variation  in  either 
eccentricity  or  inclination.  The  location  of  these 
inclinations  through  numerical  means  was  the  desired  goal. 

To  provide  a  bound  on  the  problem,  an  analytic  method 
was  outlined  (e.g.,  see  Breakwell  ref.  5).  This  resulted  in 
eleven  critical  inclinations  based  on  the  simplifications  of 
averaging  out  the  mean  anomaly  terms  and  averaging  the 
affects  of  the  Sun. 


Table  3.  Critical  Inclinations 


63.4  deg 

116.6  deg 

46.4  deg 

106.8  deg 

73.2  deg 

133.6  deg 

56.1  deg 

111.0  deg 

69.0  deg 

123.9  deg 

90.0  deg 

However,  this  method  didn't  bring  out  the  dependence  of  the 
critical  inclination  on  initial  eccentricity  and  periapsis 
height . 

The  numerical  approach  involved  the  integration  of 
orbits  over  a  range  of  initial  eccentricities  and 
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inclinations.  The  standard  deviations  in  eccentricity  (SDE) 
and  inclination  (SDI) ,  based  upon  a  polynomial  curve  fit  to 
the  data,  were  employed  as  the  search  parameters.  The  local 
maximums  in  the  surface  plots  of  SDE  and  SDI  vs  initial 
conditions  indicated  the  resonances  and  thus  the  critical 
inclinations  (e.g.,  see  figures  6-25). 

The  SDI  plots  proved  to  be  relatively  smooth  and  no 
predominate  features  were  discerned.  Therefore,  critical 
inclinations  for  large  variations  in  inclination  were  not 
determined. 

However,  the  SDE  plots  demonstrated  numerous  resonant 
conditions.  The  curves  for  critical  inclinations  were 
identified  and  roughly  corresponded  to  the  analytic  results 
(e.g.,  see  table  1). 

Experimenting  with  higher  order  polynomial  curve  fits, 
in  calculating  the  standard  deviations,  resulted  in  a 
smooching  of  the  surface  plots  as  expected.  The  more 
accurate  curve  fits  reduced  the  magnitude  of  the  deviations 
and  brought  out  the  predominate  resonance  effects  (e.g.,  see 
figure  42).  The  equations  for  these  curves  indicated  the 
critical  inclinations  were  a  function  of  the  initial 
conditions . 

A  second  effect  of  the  higher  degree  polynomial  data 
fits  was  a  shifting  of  che  critical  inclination  curves. 
Increasing  the  degree  rf  the  polynomial  shifted  the  critical 
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curves  toward  decreasing  eccentricity  and  away  from  a 
dividing  line  of  53.3  degrees  initial  inclination.  This 
dividing  line  also  corresponded  to  the  approximate 
inclination  of  the  predominate  resonance  effects  (e.g.,  see 
figure  42)  and  the  analytic  critical  inclination  of  56.1 
degrees . 

Based  on  these  results,  the  numeric  approach  utilizing 
the  standard  deviation  as  a  search  parameter  provided  a 
means  of  identifying  the  critical  inclinations  for  large 
variations  in  eccentricity.  Furthermore,  the  use  of  higher 
degree  polynomial  data  fits  in  determining  standard 
deviations  assisted  in  isolating  the  predominate  resonance 
effects . 

A  further  study  could  explore  the  shifting  of  the 
critical  inclination  curves  for  higher  degree  polynomial 
data  fits.  Non-polynomial  data  fits  may  also  be  explored. 
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(69) 


vm,n 


vm.  -a 
A0 


The  Pochhammer  symbol  (a)n  has  the  following  properties: 

U)0  =  1 

(1) a  =  n\  (70) 

(a)  n  =  a  (a- t-i)  (a+2)  ••••  (a+n-1)  n  =  1,2,3,... 

The  Pochhammer  symbol  and  the  binomial  coefficient  are 
related  by 

(a)  „  =  (-1)  n  n !  |  (71) 
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Derivation  of  Integral  Terms  in  the  Disturbing  Functions 
Derivation  of  Eqn.  (11) : 


Recognizing  this  as  a  binomial  expansion 

X03'°  =  (1  -  e2)  ^ 


(72) 


(73) 


Derivation  of  Eqn.  (11): 

Using  the  identity 

sin2f  =  —^r  [exp  (j 2f)  -  exp  ( -j2f )  ] 


(74) 
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— fcn  l  — )  sin2f  dM  - 
271  Jo  \  a  I 

=  — i-  f2"  (— )  3  — — r  [exp  (j2f)  -  exp  ( -j2f)  ]  dM 
2k  Jo  \  aj  2j 

=  f2K  (Jl\  3  exp  (j 2f)  dM  -  —  f2n  (—)  3  exp  (- j2f )  dM 

2j  [  2k  Jo  \af  2n  Jo  \  a  /  , 

=  -JJ  (*o3'2  -  *3"2)  =  0 

(75) 


Derivation  of  Eqn.  (11): 

Using  the  identity 

cos  2  f  =  [exp  ( j2f )  +  exp  ( -j2f )  ]  (7  6) 


then 

— —  [2*  (  — )  3  cos2f  dM  = 

2k  Jo  \  a) 

=  _JL  f2n  IS)  3  A  [  exp  (j2f)  +  exp  ( -  j  2f)  ]  dM 
2  k  Jo  \  a  j  2 

=  JL[_L  r*  (_£]  "  exp  ( j2f )  dM  +  —  /"2n  (  — )  3  exp  (~j2f)  dM 
2  [  2 K  Jo  \  a  J  2k  Jo  U  Z  ; 

-  j  (^o3'2  +  ^o3,2)  =  *o3'2 

(77) 


Xr 


■3 , 2 


=  0 


(78) 
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Derivation  of  Eqn.  (41)  : 


— f2K  I  — )  cos2f  dM  = 

2k  Jo  \aj 

=  -i-  f2n  (  —  j2  —  [exp(j2f)  +  exp(-j'2f)  ]  dw 
2ti  Jo  \  a  /  2 

=  r2K  /-£\2  exp  (j2f)  dM  +  -A-  f271  (— f  exp  (~j2f)  dM 

2  [  2ti  Jo  \  a)  2k  Jo  \  a  / 

=  (x02-2  +  X02-2)  =  X02'2 

(79) 


(80) 


Derivation  of  Eqn.  (41)  : 

1  <•  2  n  /  7-  \  2 

-  /  —  sin2f  dM  = 

2k  Jo  \  a  I 

=  — f2n  i  —  Y  — ~~r  ( exp  (j2f)  -  exp(-j‘2f)]  dM 
2k  Jo  \  a  I  2j 

=  i[h  f  (if exp  ^ "  -h  f  (i f exp  (-^f)  ^  j 

=  (a:2-2  -  x2--2)  =  o 

(81) 
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Derivation  of  Eqn.  (41) : 


,2,0 

v0 


£  r  (ir  ■  * 

(-T)n  (_1)  n  e2n 


=  (1)  (1)  £ 


72  =  0 


(1) 


n ! 


1  +  i-e2 
2 


(82) 
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Figure  43.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  5500  Km, 
Linear  Data  Fit 
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Figure  44.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  5500  Km, 
Quadratic  Data  Fit 


Figure  45.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  GuOu  Km, 
Linear  Data  Fit 
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Figure  46.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  6000  Km, 


Figure  47.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  6500  Km, 
Linear  Data  Fit 
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Figure  48.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  6500  Km, 
Quadratic  Data  Fit 


Figure  49.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  7000  Km, 
Linear  Data  Fit 


Figure  50.  Standard  Deviation  in  Eccentricity  vs 
Inclination  and  Eccentricity:  Periapse  Radius  =  7000  Km, 
Quadratic  Data  Fit 


Appendix  C.  SDI  Surface  Plots 


Figure  51.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5500  Km,  Linear  Data 
Fit 
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Figure  52.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  5500  Km,  Quadratic  Data 
Fit 
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Figure  53.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  6000  Km,  Linear  Data 
Fit 


Figure  54.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  6000  Km,  Quadratic  Data 
Fit 


90 


Figure  56.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  6500  Km,  Quadratic  Data 
Fit 
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Figure  58.  Standard  Deviation  in  Inclination  vs  Inclination 
and  Eccentricity:  Periapse  Radius  =  7000  Km,  Quadratic  Data 
Fit 
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Appendix  D.  Critical  Inclination  Plots 
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